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The global phase behavior of the lattice restricted primitive model with nearest neighbor exclusion 
has been studied by grand canonical Monte Carlo simulations. The phase diagram is dominated 
by a fluid (or charge-disordered solid) to charge-ordered solid transition that terminates at the 
maximum density, p'^ax = "s/S and reduced temperature T* ~ 0.29. At that point, there is a first- 
order phase transition between two phases of the same density, one charge-ordered and the other 
charge-disordered. The liquid-vapor transition for the model is metastable, lying entirely within the 
fluid-solid phase envelope. 



INTRODUCTION 



Recent theoretical and simulation studies of ionic sys- 
tems have improved our understanding of their structure 
and thermodynamic9i'2iMiSiSiL£. One of the most suc- 
cessful and simplest ionic models is the restricted primi- 
tive model (RPM), in which the ions are viewed as equi- 
sized hard spheres carrying positive and negative charges 
of the same magnitude. The RPM exhibits vapor-liquid 
phase separation and the corresponding critical point was 
confirmed to belong to the three-dimensional Ising uni- 
versality classSiili. 

The discretized version of the RPM, the lattice re- 
stricted primitive model (LRPM), has also been exten- 
sively studied by both simulatiorkiiii24iiiiii^ii& and the- 
oretical approache3iiiiii2*i2iS2i2ii2a. In the LRPM, po- 
sitions of the positive and negative ions of diameter a 
are restricted to the sites of an underlying simple cubic 
lattice of spacing l; the parameter C, = u/l specifies how 
closely the lattice system approaches the continuum be- 
havior. In addition to the obvious computational advan- 
tagesi^, since the interactions between all lattice sites are 
pre-computed, the lattice model presents some unusual 
characteristics generally absent from non-ionic fluidsS^. 
For C = Ij the most striking feature is an order-disorder 



transitio: 



11.12 



, which is not present in the continuous ver- 
sion of the RPM. There is no gas-liquid transition and 
the coexistence is between a low-density disordered phase 
and an antiferromagnetically ordered high-density phase. 
The transition is continuous (Neel-type line) above and 
first-order below a tricritical point. However, for fine dis- 
cretized lattices with C ^ 3 the vapor-liquid coexistence 
is recovered and the critical point and coexistence curves 
converge to the values found in the continuum model for 
increasing values of (^iLiSiS,*. 
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Several investigations have been performed of the C = 1 
LRPM with additional short-range attractionsi^iii or 
nearest-neighbor (nn) repulsionsiSiiSiiS. These models 
present a rich phase behavior as the nn strength is varied. 
For weak or vanishing nn interactions only order-disorder 
transitions and a tricritical point were found, while for 
increasing nn strength different scenarios could be possi- 
ble, depending on the magnitude of the nn interaction. 
For short-range attractions both tricritical and gas-liquid 
critical points can become stable^^, while for nn repulsion 
the continuum-space behavior is recoveredi^. 

The present paper extends the previous studjsiSi, of 
the LRPM model with variable nearest neighbor repul- 
sion to the limit of infinite repulsion. The model is thus 
the lattice restricted primitive model with nearest neigh- 
bor exclusion (LRPM-nn). This is equivalent to a lat- 
tice restricted primitive model of discretization param- 
eter Q = \/2. Prior theoretical studiesiSiiS with first, 
second and third nn exclusion indicate close similarity 
between the LRPM-nn phase diagram and the results 
for an off-lattice ionic system at high densities. Of par- 
ticular interest to the present work are possible connec- 
tions between order-disorder transitions and lower den- 
sity vapor-liquid transitions and high density transitions 
between charge-ordered and charge-disordered phasesi^. 
The present paper is organized as follows. The model 
and computational details are given in Sec. m Results 
are discussed in Sec. IIIII We close in Sec. II VI with sum- 
mary and conclusions. 



II. MODEL AND SIMULATION METHODS 

We consider a system of 2N charges, half of them car- 
rying charge +q and half charge — g, on a simple cu- 
bic lattice. We enforce nearest neighbor exclusion on 
the lattice and set the charge diameter as the unit of 
length, a — 1. The lattice spacing is then I — l/\/2, such 
that the lattice discretization parameter^^ is defined as 
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C = a/l = The charges interact through the (con- 
tinuous space) Coulomb potential 



^ Dr, 



(1) 



where D is the dielectric constant of the structureless 
solvent in which the charges are immersed. While some 
studies of the phase behavior of the true lattice Coulomb 
potential are available— , most prior studies of lattice 
RPM models have been done using the potential of equa- 
tion 

Reduced quantities are defined as follows: 



Eo ' 



P = 



V 



rr* U 



(2) 



where a is the ion diameter (taken as the unit of length), 
V is the volume of the system, U the energy per ion pair 
and Eq = /Da is the magnitude of Coulomb energy 
between two ions at close contact. The reduced chemical 
potential, /i* , is defined so that at the limit of high tem- 
peratures and low densities, ji* —* 2T* \nNa^ /V, where 
the factor 2 comes from the presence of two ions per min- 
imal neutral "molecule" inserted or deleted in the simu- 
lations. With this choice of reduced units, the maximum 
density of the system at close packing is pJnax = v^- 

The reduced box length, L* = L/a is a non-integer 
quantity, as the lattice spacing is I = 1/V2; it is more 
convenient to use the lattice spacing as the reducing 
length in this case, so we define 



/ a 



(3) 



The LRPM-nn model studied in the present work is 
equivalent to the limiting case of infinite repulsive cou- 
pling to nearest neighbor sites, J ^ oo, in the model 
studied in Ref.^^. However, Refi^^ used an effective 
charge diameter of u = -^72, so that reduced densities 
and temperatures are higher in the present study by fac- 
tors of 2\/2 and \/2, respectively. The unit conventions 
used here facilitate comparisons with previous studies of 
the LRPM model^^ and the continuous RPM-?MI. 

Electrostatic interactions were computed using Ewald 
summation with conducting boundary conditions at in- 
finity, 518 Fourier-space wave vectors and real-space 
damping parameter n — 5. Interactions were pre- 
computed at the beginning of the simulation runs for all 
possible pairs of lattice sites and stored in an array for 
computational efficiency. 

We used grand canonical Monte Carlo (GCMC) simu- 
lations in cubic boxes subject to full periodic boundary 
conditions. Two types of moves were utilized, namely 
pair additions and removals and swapping of oppositely 
charged ions to enhance sampling of order-disorder tran- 
sitions. To enhance acceptance of the insertion and re- 
moval steps we used distance-biased sampling, following 
Ref.^^. Particle swaps constituted up to 60 % of at- 
tempted moves, depending on temperature and density. 




FIG. 1: (color online) Normalized density distribution, F{p*), 
at T* = 0.15 in simulation boxes of size — 12 (continuous 
line) and = 8 (dashed line). 



Multihistogram reweighting2422i2£ techniques were 
used to analyze the simulation data. For the critical re- 
gion we used mixed-field finite size scaling analysis pro- 
posed by Bruce and WildingSi, which accounts for the 
lack of symmetry between coexisting phases in fluids. We 
did not attempt to incorporate corrections for pressure 
mixing in the scaling fields, as any such effects are ex- 
pected to be small^S. The Tsypin and Blote'^'^ limiting 
distribution for the three-dimensional Ising model was 
used for obtaining the critical parameters. Typical runs 
involved 10* Monte Carlo steps for equilibration and 10^ 
steps for production. Such runs required approximately 
10 CPU hours on 3 GHz Pentium 4 processors. Longer 
runs were performed near the vapor-liquid critical point 
and at high densities at which the acceptance ratio for 
insertions and removals was lower. Statistical uncertain- 
ties were obtained from multiple independent runs with 
different pseudorandom sequences. The random number 
generator "ran2" of Ref.'^ was used for the calculations. 

For the charge-disordered solid to charge-ordered solid 
part of the phase diagram (or "fluid" -solid), we were able 
to observe direct transitions between the ordered and dis- 
ordered phases, as shown in Fig. This established the 
relative free energies of the two phases and eliminated 
the need for thermodynamic integration to a reference 
state of known free energy. As seen in Fig.^ for smaller 
boxes the probability of intermediate densities is greatly 
enhanced and the system readily passes between the two 
phases. For simulation boxes much greater than — 12, 
sampling is restricted to the phase from which the simula- 
tion is started; it is not possible to establish a reversible 
path from the charge-disordered solid (or fluid) to the 
charge-ordered solid phases using grand canonical simu- 
lations without umbrella sampling. 
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FIG. 2: (color online) Phase behavior of the LRPM-nn model. 
The "fluid" -solid transition is indicated by crosses (L^ — f 2) 
and open circles (L^ = 8). Filled circles mark the vapor- 
liquid equilibrium phase boundaries obtained in a box of 
size = 24. The inset expands the region around the 
metastable vapor- liquid critical point. Statistical uncertain- 
ties are smaller than symbol size. 



III. RESULTS AND DISCUSSION 

Fig-Elshows the overall phase behavior for the LPRM- 
nn model. The "fluid" -solid transition dominates the 
phase diagram. It should be pointed out that the desig- 
nation "fluid" is not appropriate at high densities. There 
is no first-order transition between fluid and solid phases 
in the V2 lattice hard sphere model to which the present 
model reduces at the limit of high temperatures. The V2 
lattice hard sphere model has a second order transition 
at a density of p* ~ 0.59^^1^, above which an fcc-ordered 
sublattice exists in the system. The lower-density phase 
for the LPRM-nn is a disordered fluid at lower tempera- 
tures but becomes increasingly solid-like at higher densi- 
ties. Near the end-point of the transition the coexisting 
phases are both solids with face-centered-cubic overall ar- 
rangement of the particles, if one ignores the ion charge. 
For this reason, we enclose the term "fluid" in quotation 
marks when referring to the phase at the low-density side 
of this transition, to acknowledge the fact that the phase 
changes continuously from a disordered fluid to a charge- 
disordered solid. 

The vapor- liquid envelope is seen in the inset to Fig. [21 
to be metastable and wholly within the "fluid" -solid 
boundary. This is surprising for a system with long-range 
(Coulombic) interactions because the lack of a stable liq- 
uid phase is usually associated with short-range attrac- 
tions^i. This behavior can be rationalized by considering 
that the presence of the underlying lattice greatly stabi- 




FIG. 3: (color online) . Order parameter distribution, F(x), 
for a system of = 24 at conditions corresponding to the 
critical point listed in Table I. Points are simulation data and 
the line represents the analytical approximation of Refi— for 
the three-dimensional Ising universality class limiting distri- 
bution. Simulation statistical uncertainties are comparable to 
symbol size. 



lizes the solid, thus shifting the fluid-solid transition to 
higher temperatures. By contrast, for systems with very 
short range interactions, it is the destabilization of the 
liquid phase that leads to the disappearance of the vapor- 
liquid transition. The phase diagram for the LRPM-nn 
model is similar to that of the LRPM with ^ = 2, ob- 
tained with large statistical uncertainties in Refi^ and 
studied theoretically inSS. 

There is some system size dependence observed for 
the "fluid" -solid part of the phase diagram, especially 
at lower temperatures, as seen from the difference of the 
apparent phase boundaries. The smallest of the system 
sizes shown in Fig.|21(it = 8) can accommodate only 256 
ions at full packing; it was included only for comparison 
purposes. The larger system (L^ = 12) can accommo- 
date 864 ions at close packing, but is still small enough 
to allow for efficient direct sampling of the order-disorder 
transition, as seen in Fig. ^ 

The liquid-vapor critical point and phase coexistence 
envelope were obtained for system sizes ranging from 
L'l' = 14 to — 24. Because of the lower densities rel- 
ative to the "fluid" -solid transition, larger systems were 
required. Results for the critical parameters are shown 
in Table I. The apparent (system-size dependent) crit- 
ical temperature, T*, chemical potential, /i* and held 
mixing parameter s^i^ were obtained by minimizing de- 
viations between the universal order parameter distribu- 
tionS and the observed distributions. A typical opti- 
mized order parameter distribution is shown in Fig. |21 , 
in which the abscissa is the mixed-field order parameter, 
X = N{U* — s^ix). The estimates for the critical temper- 
ature do not vary significantly with system size, while the 
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TABLE L Vapor-liquid critical parameters. Numbers in 
parentheses indicate statistical uncertainties in units of the 
last figure shown. 





T* 




^mix Pc 


15 


0.0623(1) 


-1.5015(1) 


-0.643(1) 0.132(1) 


18 


0.0624(1) 


-1.5021(1) 


-0.621(1) 0.155(2) 


22 


0.0622(1) 


-1.5016(2) 


-0.630(5) 0.150(4) 


24 


0.0622(1) 


-1.5041(1) 


-0.629(3) 0.148(3) 




critical density first increases and then becomes slightly 
lower with system size. 

These critical parameters can be compared with those 
for the related model of Ref.^^ with partial exclusion of 
nearest neighbor sites, J ~ 0.3, after taking into account 
the different reducing parameters. The ,/ = 0.3 model 
has critical parameters (for it = 15) T* = 0.0612(1), 
p* ~ 0.14(2), very similar to the LPRM-nn present model 
which corresponds to J — > oo. However, the high-density 
(solid) phases of the two models are completely different, 
as discussed in the paragraphs that follow. 

The critical temperature for the present model is 
slightly higher and the critical density significantly higher 
than the continuum RPM, for which T* = 0.04933(5), 
p* = 0.075(l)i2'. Vapor-liquid critical parameters for 
models with ^ = 3, 4 and 5 have been obtained ir*i^; 
the critical temperatures and densities were found to be 
higher in coarser lattices, a trend consistent with the re- 
sults of the present study. 

At relatively high temperatures, Coulombic interac- 
tions become less important than excluded volume in de- 
termining the equation of state for the LPRM-nn model. 
At the limit of very high temperatures, we have already 
mentioned that the model is equivalent to a \/2 lattice 
hard sphere model. At higher densities an fcc-ordered 
solid appears and at close packing every sphere occupies 
an ordered position. In LRPM-nn, the high-temperature 
solid is substitutionally disordered with respect to charge 
type. The continuous RPM model^^'^^''^^ has a transi- 
tion from a charge-disordered to a charge-ordered solid 
phase at temperatures near T* k, 0.29. The fully oc- 
cupied LRPM charge-ordered to charge-disordered phase 
transition has been investigated previously'^^ and found 
to have a first order transition at the same temperature. 
The phase behavior of Fig. ^ has the same transition at 
a similar temperature range. 

The structure of the ion-ordered solid is shown in 
Fig. El The structure is P4/mmm (tetragonal), identi- 
cal to the "fee" -ordered structure observed at high den- 
sities for the continuum RPM^SiSLS. This structure has 
not yet been observed experimentally in systems of op- 
positely charged colloids^. 

In Ref4S, the character of the transition between or- 
dered and disordered solid phases for the continuous 
RPM was investigated using constant-pressure Monte 
Carlo simulations and found to be "weakly first order." 




FIG. 4: (color online) Structure of the solid, = 12. Differ- 
ent shades (colors) represent oppositely charged ions. 
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FIG. 5: (color online) Density (top) and energy (bottom) 
distributions, _F(p*) and F{U*) at T* = 0.28, /i* = 0.8. Solid 
line: — 12, started from disordered solid; dashed line: 
I/t = 12 started from ordered solid; dash-dotted line: — 8. 



For the LRPM-nn model, the densities of the coexist- 
ing phases (ordered and disordered solids) are seen in 
Fig. 121 to converge at the maximum density p^Q^, = \/2. 
Normalized probability distributions for the densities and 
energies are shown in Fig. Oat T* — 0.28 for two system 
sizes, = 12 and — xhe apparent "noise" at low 
energies (to the right of the bottom part of the graph) is 
due to the finite number of states with one, two etc vacan- 
cies in the lattice model. For the smaller box, densities 
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FIG. 6: (color online) Density difference between charged- 
ordered and charged-disordered phases for = 12 (crosses) 
and — 8 (open circles). Statistical uncertainties are com- 
parable to symbol size. Lines are linear-least-squares fits to 
the points over the range indicated. 



and energies corresponding to both phases are sampled 
in a single run. There is no hint of two peaks in the 
density distribution, but the energy distribution shows a 
clear separation of states into higher energy (less nega- 
tive, disordered) and lower energies (ordered). For the 
larger box size, the simulations get trapped in the phase 
from which they are started; even though the density dis- 
tributions of the two runs at identical chemical potentials 
overlap at p* « 1.39, there is no overlap in the energy 
distributions. 

We have computed the density difference between the 
charged-ordered and charged-disordered phases near the 
transition end-point by using energy (rather than den- 
sity) to identify the phases. In other words, referring 
to Fig. [SI we collected the density under the charge- 
ordered (more negative energy) and charge-disordered 
(less negative energy) phases. The results allow us to 
extend the coexistence envelope to higher temperatures 
for which the density differences are small. The results 
for the density difference, Ap*, as a function of tempera- 
ture are shown in Fig. El The linear-least-squares fits to 
the points are indicated as lines in Fig. El These give in 
turn the temperature at which the first order transition 
occurs with no density discontinuity, as T* = 0.295(2) 
for Lt = 8 and T* ^ 0.291(3) for = 12. At that 
limit, both coexisting phases are at the closed-packed 
density, p* = p^ax — V^- ^ fluid-solid phase transi- 
tion with no density discontinuity has been established 
for the Gaussian core model^. Experimentally, several 
metallic elements, in particular Ce, Cs, Ba and Eu, have 
melting lines of zero slope in the pressure-temperature 
pjg^jjg42j43j44^ also indicating a fluid-solid transition with 
no density discontinuity. 

It is of interest to compare our results to the field- 



theoretic study of Ciach and StelUS for the LPRM -nn and 
to simulations of Abascal et alm^ and Bresme et al?^ for 
the continuous RPM at high densities. The main differ- 
ence between our results and these prior studies is that 
we find that the charge-ordered and charge-disordered 
solids have a first-order transition with zero density dif- 
ference at the closed-packed density. We speculate that 
such a transition may also be found in the continuous 
RPM near the close-packed density, as our findings seem 
to be generally consistent with those of Ref.'*''. 



IV. CONCLUSIONS 

We have studied the phase diagram of the lattice re- 
stricted primitive model with nearest neighbor exclusion 
(LPRM-nn), using grand canonical Monte Carlo simula- 
tion and histogram reweighting techniques. The global 
phase diagram is dominated by the "fluid" -solid transi- 
tion, which starts with a large density gap between a 
dilute gas phase and the solid at low temperatures. The 
transition ends at T* = 0.291(3) as a first order transition 
between charge-ordered and charge-disordered phases of 
the same density, p* = \f2. First-order transitions be- 
tween phases of the same density in one-component sys- 
tems have been observed for several metallic elements 
and for the Gaussian core model. 

The liquid-vapor phase transition for the model was 
determined to be metastable, lying entirely within the 
"fluid" -solid phase envelope. The metastable critical 
point for the transition was obtained as — 0.0622(1), 
p* = 0.148(3), values higher than for the continuum 
RPM but consistent with previously determined trends 
for discretized lattice models. 

While the broad outline of the phase diagram is consis- 
tent with theoretical predictions^^, our results differ from 
these predictions in some important aspects. In particu- 
lar, the liquid / fcc-disordered transition is not present in 
our system. Our results are in near-quantitative agree- 
ment with calculations of ordered-disordered fee phase 
transitions for the continuum RPM'*^. However, we flnd 
that the charge-ordered and charge-disordered phases 
maintain a first-order transition even though there is no 
density difference between the coexisting phases. 
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